﻿//用于火力分配计算

#include <iostream>
#include <sstream>
#include <iomanip>
#include <vector>
#include <string>
#include <cmath>
#include <algorithm>

using namespace std;

//定义WGS84椭球体基准大地参数
const double a = 6378137.0; //长半轴
const double b = 6356752.314245; //短半轴
const double f = (a - b) / a; //扁率 1/298.257223563
const double e = sqrt(2 * f - pow(f, 2)); //第一偏心率
const double e2 = sqrt(pow(a, 2) - pow(b, 2)) / b; //第二偏心率
const double e12 = pow(e, 2) / (1 - pow(e, 2)); //第一偏心率平方
const double e22 = pow(e2, 2) / (1 - pow(e2, 2)); //第二偏心率平方
const double m0 = a * (1 - e12);  // 子午线曲率半径
const double n0 = a / sqrt(1 - e12 * sin(0) * sin(0));  // 卯酉圈曲率半径

const double M_PI = 3.14159265358979323846;

//地理坐标类
class GeoCoordinateStruct {
public:
    GeoCoordinateStruct() = default;
    //构造函数, 高度默认为0
    GeoCoordinateStruct(double longitude, double latitude, double altitude = 0) {
        m_longitude = longitude;
        m_latitude = latitude;
        m_altitude = altitude;

        //非法坐标处理
        if (!this->isValid()) {
            m_longitude = 0;
            m_latitude = 0;
            m_altitude = 0;
        }
    }

    //get和set方法
    double getLongitude() { return m_longitude; }
    double getLatitude() { return m_latitude; }
    double getAltitude() { return m_altitude; }
    void setLongitude(double longitude) { m_longitude = longitude; }
    void setLatitude(double latitude) { m_latitude = latitude; }
    void setAltitude(double altitude) { m_altitude = altitude; }

    //是否是有效坐标
    bool isValid() {
        if (m_longitude < -180 || m_longitude > 180 || m_latitude < -90 || m_latitude > 90) {
            return false;
        }
        else {
            return true;
        }
    }

    //使用Vincenty公式计算与另一坐标点之间的距离
    double calculateDistance(GeoCoordinateStruct coordinate) {
        //非法坐标处理
        if (!this->isValid() || !coordinate.isValid())
            return 0;

        GeoCoordinateStruct coordinate1(m_longitude, m_latitude);
        GeoCoordinateStruct coordinate2(coordinate.getLongitude(), coordinate.getLatitude());
        return GeoCoordinateStruct::calculateDistance(coordinate1, coordinate2);
    }

    //使用Vincenty公式计算两点之间的距离，精度0.001米
    static double calculateDistance(GeoCoordinateStruct coordinate1, GeoCoordinateStruct coordinate2) {
        //非法坐标处理
        if (!coordinate1.isValid() || !coordinate2.isValid())
            return 0;

        double lon1 = toRadians(coordinate1.getLongitude());
        double lat1 = toRadians(coordinate1.getLatitude());
        double lon2 = toRadians(coordinate2.getLongitude());
        double lat2 = toRadians(coordinate2.getLatitude());

        double L = lon2 - lon1;
        double U1 = atan((1 - f) * tan(lat1));
        double U2 = atan((1 - f) * tan(lat2));
        double sinU1 = sin(U1);
        double cosU1 = cos(U1);
        double sinU2 = sin(U2);
        double cosU2 = cos(U2);

        double lambda = L;
        double lambdaP = 2 * M_PI;
        double iterLimit = 20;
        double cosSqAlpha = 0;
        double sinSigma = 0;
        double cos2SigmaM = 0;
        double cosSigma = 0;
        double sigma = 0;
        double sinAlpha = 0;
        double C = 0;

        while (fabs(lambda - lambdaP) > 1e-12 && --iterLimit > 0) {
            double sinLambda = sin(lambda);
            double cosLambda = cos(lambda);
            sinSigma = sqrt(pow(cosU2 * sinLambda, 2) + pow(cosU1 * sinU2 - sinU1 * cosU2 * cosLambda, 2));
            if (sinSigma == 0) {
                return 0;
            }
            cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
            sigma = atan2(sinSigma, cosSigma);
            sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;
            cosSqAlpha = 1 - pow(sinAlpha, 2);
            cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
            if (isnan(cos2SigmaM)) {
                cos2SigmaM = 0;
            }
            C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha));
            lambdaP = lambda;
            lambda = L + (1 - C) * f * sinAlpha * (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * pow(cos2SigmaM, 2))));
        }

        if (iterLimit == 0) {
            return 0;
        }

        double uSq = cosSqAlpha * (pow(a, 2) - pow(b, 2)) / pow(b, 2);
        double A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)));
        double B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)));
        double deltaSigma = B * sinSigma * (cos2SigmaM + B / 4 * (cosSigma * (-1 + 2 * pow(cos2SigmaM, 2)) - B / 6 * cos2SigmaM * (-3 + 4 * pow(sinSigma, 2)) * (-3 + 4 * pow(cos2SigmaM, 2))));
        double s = b * A * (sigma - deltaSigma);

        return s;
    }

    //计算与另一坐标点之间的方位角
    double calculateBearing(GeoCoordinateStruct coordinate) {
        double lon1 = toRadians(m_longitude);
        double lat1 = toRadians(m_latitude);
        double lon2 = toRadians(coordinate.getLongitude());
        double lat2 = toRadians(coordinate.getLatitude());

        double y = sin(lon2 - lon1) * cos(lat2);
        double x = cos(lat1) * sin(lat2) - sin(lat1) * cos(lat2) * cos(lon2 - lon1);
        double bearing = atan2(y, x);

        return fmod((toDegrees(bearing) + 360), 360);
    }

    //计算两点之间的方位角
    static double calculateBearing(GeoCoordinateStruct coordinate1, GeoCoordinateStruct coordinate2) {
        double lon1 = toRadians(coordinate1.getLongitude());
        double lat1 = toRadians(coordinate1.getLatitude());
        double lon2 = toRadians(coordinate2.getLongitude());
        double lat2 = toRadians(coordinate2.getLatitude());

        double y = sin(lon2 - lon1) * cos(lat2);
        double x = cos(lat1) * sin(lat2) - sin(lat1) * cos(lat2) * cos(lon2 - lon1);
        double bearing = atan2(y, x);

    }

    //根据方位角和距离计算另一点的坐标
    GeoCoordinateStruct calculateCoordinate(double bearing, double distance) {
        //非法坐标处理
		if (!this->isValid())
			return GeoCoordinateStruct(0, 0);

		//计算另一点的坐标
		double lon1 = toRadians(m_longitude);
		double lat1 = toRadians(m_latitude);
		double bearingRadians = toRadians(bearing);
		double angularDistance = distance / m0;
		double lat2 = asin(sin(lat1) * cos(angularDistance) + cos(lat1) * sin(angularDistance) * cos(bearingRadians));
		double lon2 = lon1 + atan2(sin(bearingRadians) * sin(angularDistance) * cos(lat1), cos(angularDistance) - sin(lat1) * sin(lat2));
		lon2 = fmod((lon2 + 3 * M_PI), (2 * M_PI)) - M_PI;

		return GeoCoordinateStruct(toDegrees(lon2), toDegrees(lat2));
    }

    //根据方位角和距离计算另一点的坐标

    //转换成度分秒格式字符串，保留小数点后三位
    std::string toDMSString() {
        double longitude = m_longitude;
        double latitude = m_latitude;

        std::ostringstream oss;

        int lonDegrees = static_cast<int>(longitude);
        int lonMinutes = static_cast<int>((longitude - lonDegrees) * 60);
        double lonSeconds = (longitude - lonDegrees - lonMinutes / 60.0) * 3600;

        int latDegrees = static_cast<int>(latitude);
        int latMinutes = static_cast<int>((latitude - latDegrees) * 60);
        double latSeconds = (latitude - latDegrees - latMinutes / 60.0) * 3600;

        oss << std::fixed << std::setprecision(2);
        oss << lonDegrees << "°" << lonMinutes << "'" << lonSeconds << "\"";
        oss << " ";
        oss << latDegrees << "°" << latMinutes << "'" << latSeconds << "\"";

        return oss.str();
    }

    //角度转换为弧度
    static double toRadians(double degrees) { return degrees * M_PI / 180; }
    //弧度转换为角度
    static double toDegrees(double radians) { return radians * 180 / M_PI; }

private:
    double m_longitude; // 经度
    double m_latitude;  // 纬度
    double m_altitude;  // 高度
};

//计算散列坐标的中心点
GeoCoordinateStruct calculateCoordsCenter(std::vector<GeoCoordinateStruct> coordinates) {
	double longitude = 0;
	double latitude = 0;
    for (auto c : coordinates) {
		longitude += c.getLongitude();
		latitude += c.getLatitude();
	}
	longitude /= coordinates.size();
	latitude /= coordinates.size();

	return GeoCoordinateStruct(longitude, latitude);
}

//直线类
class LineStruct {
public:
    //默认构造函数
    LineStruct() = default;

    //构造函数,一般式构造直线:Ax+By+C=0(A、B不同时为0)
    LineStruct(double a, double b, double c) {
        //非法直线处理
        if (a == 0 && b == 0) {
            this->invalidLine();

            return;
        }

        //计算直线的A、B、C系数
        m_a = a;
        m_b = b;
        m_c = c;

        //计算直线的斜率
        m_slope = -a / b;

        //计算直线的x截距
        m_xintercept = -c / a;

        //计算直线的y截距
        m_yintercept = -c / b;
    }

    //构造函数，点斜式构造直线:y-y0=k(x-x0)(k不为0)
    LineStruct(double slope, GeoCoordinateStruct point) {
        //非法直线处理
        if (slope == 0) {
            this->invalidLine();

            return;
        }

        //计算直线的A、B、C系数
        m_a = -slope;
        m_b = 1;
        m_c = slope * point.getLongitude() - point.getLatitude();

        //计算直线的斜率
        m_slope = slope;

        //计算直线的x截距
        m_xintercept = -m_c / m_a;

        //计算直线的y截距
        m_yintercept = -m_c / m_b;
    }

    //构造函数，两点式构造直线:(y-y1)/(y2-y1)=(x-x1)/(x2-x1)(x1不等于x2)
    LineStruct(GeoCoordinateStruct point1, GeoCoordinateStruct point2) {
        //非法直线处理, 两点重合或者点不合法
        if (point1.getLongitude() == point2.getLongitude() || !point1.isValid() || !point2.isValid()) {
            this->invalidLine();

            return;
        }

        //计算直线的A、B、C系数
        m_a = point2.getLatitude() - point1.getLatitude();
        m_b = point1.getLongitude() - point2.getLongitude();
        m_c = point2.getLongitude() * point1.getLatitude() - point1.getLongitude() * point2.getLatitude();

        //计算直线的斜率
        m_slope = -m_a / m_b;

        //计算直线的x截距
        m_xintercept = -m_c / m_a;

        //计算直线的y截距
        m_yintercept = -m_c / m_b;
    }

    //构造函数，斜截式构造直线:y=kx+b(k不为0)
    LineStruct(double slope, double intercept) {
        //非法直线处理
        if (slope == 0) {
            this->invalidLine();

            return;
        }

        //计算直线的A、B、C系数
        m_a = -slope;
        m_b = 1;
        m_c = -intercept;

        //计算直线的斜率
        m_slope = slope;

        //计算直线的x截距
        m_xintercept = -m_c / m_a;

        //计算直线的y截距
        m_yintercept = -m_c / m_b;
    }

    //非法直线参数处理
    void invalidLine() {
        m_a = 0;
        m_b = 0;
        m_c = 0;

        m_slope = 0;

        m_xintercept = 0;
        m_yintercept = 0;
    }

    //get方法
    double getA() { return m_a; }
    double getB() { return m_b; }
    double getC() { return m_c; }

    double getSlope() { return m_slope; }
    double getXIntercept() { return m_xintercept; }
    double getYIntercept() { return m_yintercept; }

    //计算垂直于直线的斜率
    double calculatePerpendicularSlope() {
        //非法直线处理
        if (!this->isValid())
            return 0;

        //计算垂直于直线的斜率
        double perpendicularSlope = -1 / m_slope;

        return perpendicularSlope;
    }

    //是否是有效直线
    bool isValid() {
        if (m_a == 0 && m_b == 0) {
            return false;
        }
        else {
            return true;
        }
    }

    //计算点到直线的距离
    double calculateDistance(GeoCoordinateStruct point) {
        //非法直线处理
        if (!this->isValid())
            return 0;

        //计算点到直线的距离
        double distance = fabs(m_a * point.getLongitude() + m_b * point.getLatitude() + m_c) / sqrt(pow(m_a, 2) + pow(m_b, 2));

        return distance;
    }

    //计算点到直线的垂足坐标
    GeoCoordinateStruct calculateFoot(GeoCoordinateStruct point) {
        //非法直线处理
        if (!this->isValid())
            return GeoCoordinateStruct(0, 0);

        //计算点到直线的垂足坐标
        double x = (m_b * m_b * point.getLongitude() - m_a * m_b * point.getLatitude() - m_a * m_c) / (pow(m_a, 2) + pow(m_b, 2));
        double y = (m_a * m_a * point.getLatitude() - m_a * m_b * point.getLongitude() - m_b * m_c) / (pow(m_a, 2) + pow(m_b, 2));

        return GeoCoordinateStruct(x, y);
    }

    //判断与另一条直线是否相交
    bool isIntersect(LineStruct line) {
        //非法直线处理
        if (!this->isValid() || !line.isValid())
            return false;

        //判断两条直线是否相交
        if (this->getSlope() == line.getSlope()) {
            return false;
        }
        else {
            return true;
        }
    }

    //计算与另一条直线的交点坐标
    GeoCoordinateStruct calculateIntersection(LineStruct line) {
        //非法直线处理
        if (!this->isValid() || !line.isValid())
            return GeoCoordinateStruct(0, 0);

        //计算两条直线的交点坐标
        GeoCoordinateStruct intersection = LineStruct::calculateIntersection(*this, line);
        return intersection;
    }

    //计算两条直线的交点坐标
    static GeoCoordinateStruct calculateIntersection(LineStruct line1, LineStruct line2) {
		//非法直线处理
		if (!line1.isValid() || !line2.isValid())
			return GeoCoordinateStruct(0, 0);

        //计算两条直线的交点坐标
        double x = (line1.getB() * line2.getC() - line2.getB() * line1.getC()) / (line1.getA() * line2.getB() - line2.getA() * line1.getB());
        double y = (line1.getA() * line2.getC() - line2.getA() * line1.getC()) / (line1.getB() * line2.getA() - line2.getB() * line1.getA());

        return GeoCoordinateStruct(x, y);
	}

    //计算过某点的平行线
    LineStruct calculateParallelLine(GeoCoordinateStruct point) {
        //非法直线处理
        if (!this->isValid())
            return LineStruct(0, 0, 0);

        //计算过某点的平行线
        double slope = this->getSlope();
        double intercept = point.getLatitude() - slope * point.getLongitude();

        return LineStruct(slope, intercept);
    }

    //计算过某点的垂直线
    LineStruct calculatePerpendicularLine(GeoCoordinateStruct point) {
        //非法直线处理
        if (!this->isValid())
            return LineStruct(0, 0, 0);

        //计算过某点的垂直线
        double slope = this->calculatePerpendicularSlope();
        double intercept = point.getLatitude() - slope * point.getLongitude();

        return LineStruct(slope, intercept);
    }

    //输出直线信息
    std::string toString() {
        std::stringstream ss;
        ss << "A: " << m_a << std::endl;
        ss << "B: " << m_b << std::endl;
        ss << "C: " << m_c << std::endl;

        //输出直线的斜率
        double slope = this->getSlope();
        ss << "斜率: " << slope << std::endl;

        //输出直线的x截距
        double xintercept = this->getXIntercept();
        ss << "x截距: " << xintercept << std::endl;

        //输出直线的y截距
        double yintercept = this->getXIntercept();
        ss << "y截距: " << yintercept << std::endl;

        return ss.str();
    }

private:
    double m_a;	//直线的A系数
    double m_b;	//直线的B系数
    double m_c;	//直线的C系数

    double m_slope;			//直线的斜率
    double m_xintercept;	//直线的x截距
    double m_yintercept;	//直线的y截距
};

//矩形类
class RectangleStruct {
public:
    //默认构造函数
    RectangleStruct() = default;

    //构造函数,根据中心点坐标、方位角、正面和纵深构造矩形
    //方位角为正北方向为0度，顺时针方向为正，逆时针方向为负
    //方位角对应的边为正面，方位角加180度对应的边为纵深
    RectangleStruct(GeoCoordinateStruct center, double bearing, double front, double depth) {
        //非法矩形处理
		if (!center.isValid() || front <= 0 || depth <= 0) {
			m_isValid = false;
			return;
		}

		//计算矩形的中心点坐标
		m_center = center;

		//计算矩形的方位角
		m_bearing = bearing;

		//计算矩形的正面和纵深
		m_front = front;
		m_depth = depth;

        //计算顶部和底部的中心点坐标
        GeoCoordinateStruct topCenter = m_center.calculateCoordinate(m_bearing, m_front / 2);
        GeoCoordinateStruct bottomCenter = m_center.calculateCoordinate(m_bearing + 180, m_front / 2);
#if 0   //这样计算在界面显示有问题，显示为平行四边形
        //计算矩形的四个顶点坐标
        m_topLeft = topCenter.calculateCoordinate(m_bearing - 90, m_depth / 2);
        m_topRight = topCenter.calculateCoordinate(m_bearing + 90, m_depth / 2);
        m_bottomLeft = bottomCenter.calculateCoordinate(m_bearing - 90, m_depth / 2);
        m_bottomRight = bottomCenter.calculateCoordinate(m_bearing + 90, m_depth / 2);
#endif

        //顶部和底部的中心点坐标连为直线，并计算垂直线
        LineStruct line(topCenter, bottomCenter);
        LineStruct perpendicularLine = line.calculatePerpendicularLine(center);

        //计算四条边
        LineStruct leftLine = line.calculateParallelLine(center.calculateCoordinate(m_bearing - 90, m_depth / 2));
        LineStruct rightLine = line.calculateParallelLine(center.calculateCoordinate(m_bearing + 90, m_depth / 2));
        LineStruct topLine = perpendicularLine.calculateParallelLine(center.calculateCoordinate(m_bearing, m_front / 2));
        LineStruct bottomLine = perpendicularLine.calculateParallelLine(center.calculateCoordinate(m_bearing + 180, m_front / 2));

        //计算四个顶点
        m_topLeft = leftLine.calculateIntersection(topLine);
        m_topRight = rightLine.calculateIntersection(topLine);
        m_bottomLeft = leftLine.calculateIntersection(bottomLine);
        m_bottomRight = rightLine.calculateIntersection(bottomLine);
    }

    //构造函数,根据4个点坐标构造矩形，根据参考点确定矩形的方位角、正面和纵深
    RectangleStruct(std::vector<GeoCoordinateStruct> points, GeoCoordinateStruct reference) {
        //非法参数处理
        if (points.size() != 4 || !reference.isValid()) {
            m_isValid = false;
        }

        //计算4个点的中心点坐标
        m_center = calculateCoordsCenter(points);

        //计算矩形的方位角
        m_bearing = reference.calculateBearing(m_center);   

        //找到矩形的左上顶点
        int index = 0;
        for (int i = 0; i < points.size(); i++) {
            double bearing = m_center.calculateBearing(points[i]);
			if (bearing < m_bearing) {
				index = i;
				break;
			}
        }

        //根据矩形的方位找到左上、右上、右下、左下4个点的坐标
        m_topLeft = points[index];
		m_topRight = points[(index + 1) % points.size()];
		m_bottomRight = points[(index + 2) % points.size()];
		m_bottomLeft = points[(index + 3) % points.size()];

		//计算矩形的正面和纵深
		m_front = GeoCoordinateStruct::calculateDistance(m_topLeft, m_topRight);
		m_depth = GeoCoordinateStruct::calculateDistance(m_topLeft, m_bottomLeft);
    }

    //构造函数,根据4条直线坐标构造矩形，根据参考点确定矩形的方位角、正面和纵深
    RectangleStruct(std::vector<LineStruct> lines, GeoCoordinateStruct reference) {
        //非法参数处理
        if (lines.size() != 4 || !reference.isValid()) {
			m_isValid = false;
		}

		//按顺序计算4条直线的交点坐标
        std::vector<GeoCoordinateStruct> points;
        for (int i = 0; i < lines.size(); i++) {
            GeoCoordinateStruct point = lines[i].calculateIntersection(lines[(i + 1) % lines.size()]);
			points.push_back(point);
        }

        //计算4个点的中心点坐标
        m_center = calculateCoordsCenter(points);

        //计算矩形的方位角
        m_bearing = reference.calculateBearing(m_center);

        //找到矩形的左上顶点
        int index = 0;
        for (int i = 0; i < points.size(); i++) {
	        double bearing = m_center.calculateBearing(points[i]);
	        
            if (bearing < m_bearing) {
                index = i;
                break;
            }
        }

        //根据矩形的方位找到左上、右上、右下、左下4个点的坐标
        m_topLeft = points[index];
        m_topRight = points[(index + 1) % points.size()];
        m_bottomRight = points[(index + 2) % points.size()];
        m_bottomLeft = points[(index + 3) % points.size()];

		//计算矩形的正面和纵深
		m_front = GeoCoordinateStruct::calculateDistance(m_topLeft, m_topRight);
		m_depth = GeoCoordinateStruct::calculateDistance(m_topLeft, m_bottomLeft);
    }

    //判断矩形是否有效
    bool isValid() {
		return m_isValid;
	}

    //get方法
    GeoCoordinateStruct getCenter() { return m_center; }
    double getBearing() { return m_bearing; }
    double getFront() { return m_front; }
    double getDepth() { return m_depth; }
    GeoCoordinateStruct getTopLeft() { return m_topLeft; }
    GeoCoordinateStruct getTopRight() { return m_topRight; }
    GeoCoordinateStruct getBottomRight() { return m_bottomRight; }
    GeoCoordinateStruct getBottomLeft() { return m_bottomLeft; }

    //获取矩形顶点的最大最小纬度和经度
    double getMaxLatitude() {
        double maxLatitude = m_topLeft.getLatitude();
        if (m_topRight.getLatitude() > maxLatitude) {
            maxLatitude = m_topRight.getLatitude();
        }
        if (m_bottomRight.getLatitude() > maxLatitude) {
            maxLatitude = m_bottomRight.getLatitude();
        }
        if (m_bottomLeft.getLatitude() > maxLatitude) {
            maxLatitude = m_bottomLeft.getLatitude();
        }
        return maxLatitude;
    }

    double getMinLatitude() {
        double minLatitude = m_topLeft.getLatitude();
        if (m_topRight.getLatitude() < minLatitude) {
            minLatitude = m_topRight.getLatitude();
        }
        if (m_bottomRight.getLatitude() < minLatitude) {
            minLatitude = m_bottomRight.getLatitude();
        }
        if (m_bottomLeft.getLatitude() < minLatitude) {
            minLatitude = m_bottomLeft.getLatitude();
        }
        return minLatitude;
    }

    double getMaxLongitude() {
        double maxLongitude = m_topLeft.getLongitude();
        if (m_topRight.getLongitude() > maxLongitude) {
            maxLongitude = m_topRight.getLongitude();
        }
        if (m_bottomRight.getLongitude() > maxLongitude) {
            maxLongitude = m_bottomRight.getLongitude();
        }
        if (m_bottomLeft.getLongitude() > maxLongitude) {
            maxLongitude = m_bottomLeft.getLongitude();
        }
        return maxLongitude;
    }

    double getMinLongitude() {
        double minLongitude = m_topLeft.getLongitude();
        if (m_topRight.getLongitude() < minLongitude) {
            minLongitude = m_topRight.getLongitude();
        }
        if (m_bottomRight.getLongitude() < minLongitude) {
            minLongitude = m_bottomRight.getLongitude();
        }
        if (m_bottomLeft.getLongitude() < minLongitude) {
            minLongitude = m_bottomLeft.getLongitude();
        }
        return minLongitude;
    }

    //计算矩形的中心点
    GeoCoordinateStruct calculateCenter() {
        return m_center;
    }

    /**
     * 正面和纵深分别按照m和n等分后返回所有的矩形
     * @param m 正面方向分割数
     * @param n 纵深方向分割数
     * @return 矩形向量
     */
    std::vector<RectangleStruct> split(int m, int n) {
        std::vector<RectangleStruct> rectangles;

        //取正面两条边的截距，计算正面分割步长
        LineStruct topLine(m_topLeft, m_topRight);
        LineStruct bottomLine(m_bottomLeft, m_bottomRight);

        double topIntercept = topLine.getYIntercept();
        double bottomIntercept = bottomLine.getYIntercept();
        double frontStep = (topIntercept - bottomIntercept) / m;

        //计算出正面方向的分割线
        std::vector<LineStruct> frontLines;
        for (int i = 0; i <= m; i++) {
			double intercept = bottomIntercept + i * frontStep;
            LineStruct line(bottomLine.getSlope(), intercept);
			frontLines.push_back(line);
		}

        //取纵深两条边的截距，计算纵深分割步长
		LineStruct leftLine(m_topLeft, m_bottomLeft);
		LineStruct rightLine(m_topRight, m_bottomRight);

		double leftIntercept = leftLine.getYIntercept();
		double rightIntercept = rightLine.getYIntercept();
		double depthStep = (rightIntercept - leftIntercept) / n;

		//计算出纵深方向的分割线
		std::vector<LineStruct> depthLines;
        for (int i = 0; i <= n; i++) {
            double intercept = leftIntercept + i * depthStep;
            LineStruct line(leftLine.getSlope(), intercept);
			depthLines.push_back(line);
        }

        //计算出所有的矩形
        for (int i = 0; i < frontLines.size() - 1; i++) {
			for (int j = 0; j < depthLines.size() - 1; j++) {
                std::vector<LineStruct> lines;
				lines.push_back(frontLines[i]);
				lines.push_back(depthLines[j]);
				lines.push_back(frontLines[i + 1]);
				lines.push_back(depthLines[j + 1]);
				RectangleStruct rectangle(lines, m_center);
                if (rectangle.isValid())
                    rectangles.push_back(rectangle);
			}
		}

        return rectangles;
    }


    //判断一个点是否在矩形内
    bool contains(GeoCoordinateStruct point) {
        //判断点是否合法
        if (!point.isValid()) {
			return false;
		}

        //判断点是否在矩形的四条边上
        LineStruct line1(m_topLeft, m_topRight);
        LineStruct line2(m_topRight, m_bottomRight);
        LineStruct line3(m_bottomRight, m_bottomLeft);
        LineStruct line4(m_bottomLeft, m_topLeft);

        double distance1 = GeoCoordinateStruct::calculateDistance(point, line1.calculateIntersection(LineStruct(m_topLeft, m_bottomLeft)));
        double distance2 = GeoCoordinateStruct::calculateDistance(point, line2.calculateIntersection(LineStruct(m_topRight, m_bottomRight)));
        double distance3 = GeoCoordinateStruct::calculateDistance(point, line3.calculateIntersection(LineStruct(m_bottomLeft, m_bottomRight)));
        double distance4 = GeoCoordinateStruct::calculateDistance(point, line4.calculateIntersection(LineStruct(m_topLeft, m_topRight)));

        //浮点数比较，如果两个浮点数的差值小于0.1则认为相等
        if (fabs(distance1) < 0.1 || fabs(distance2) < 0.1 || fabs(distance3) < 0.1 || fabs(distance4) < 0.1) {
            return true;
        }

        //判断点是否在矩形内部
        if (point.getLongitude() > m_topLeft.getLongitude() && point.getLongitude() < m_topRight.getLongitude() && point.getLatitude() > m_bottomLeft.getLatitude() && point.getLatitude() < m_topLeft.getLatitude()) {
            return true;
        }

        return false;
    }

    //输出矩形信息
    std::string toString() {
        std::stringstream ss;
        ss << "中心点坐标: " << m_center.toDMSString() << std::endl;
        ss << "方位角: " << m_bearing << std::endl;
        ss << "正面: " << m_front << std::endl;
        ss << "纵深: " << m_depth << std::endl;
        ss << "顶点坐标: " << std::endl;
        ss << "左上: " << m_topLeft.toDMSString() << std::endl;
        ss << "右上: " << m_topRight.toDMSString() << std::endl;
        ss << "右下: " << m_bottomRight.toDMSString() << std::endl;
        ss << "左下: " << m_bottomLeft.toDMSString() << std::endl;
        return ss.str();
    }

private:
    GeoCoordinateStruct m_center;   //中心点
    double m_bearing;   //方位角

    //顶点坐标
    GeoCoordinateStruct m_topLeft;
    GeoCoordinateStruct m_topRight;
    GeoCoordinateStruct m_bottomRight;
    GeoCoordinateStruct m_bottomLeft;

    //正面和纵深
    double m_front;
    double m_depth;

    //矩形是否有效
    bool m_isValid = true;
};

//定义多边形类
class PolygonStruct {
public:
    /**
     * 构造函数
     * @param points 多边形的点
     */
    PolygonStruct(std::vector<GeoCoordinateStruct> points) {
        m_points = points;
    }

    /**
     * 计算多边形中心点
     * @return 中心点
     */
    GeoCoordinateStruct calculateCenter() {
        double x = 0.0;
        double y = 0.0;
        double z = 0.0;
        for (auto point : m_points) {
            double lat = point.getLatitude() * M_PI / 180;
            double lon = point.getLongitude() * M_PI / 180;
            x += cos(lat) * cos(lon);
            y += cos(lat) * sin(lon);
            z += sin(lat);
        }
        double total = m_points.size();
        x /= total;
        y /= total;
        z /= total;
        double centerLon = atan2(y, x);
        double centerSqrt = sqrt(x * x + y * y);
        double centerLat = atan2(z, centerSqrt);

        return GeoCoordinateStruct(centerLon * 180 / M_PI, centerLat * 180 / M_PI);
    }

    /**
     * 判断某一坐标在不在这个区域内
     * @param point 待判断的坐标
     * @return true表示在区域内，false表示不在区域内
     */
    bool contains(GeoCoordinateStruct point) {
        int crossings = 0;
        for (int i = 0; i < m_points.size(); i++) {
            GeoCoordinateStruct a = m_points[i];
            GeoCoordinateStruct b = m_points[(i + 1) % m_points.size()];
            if (((a.getLatitude() <= point.getLatitude()) && (point.getLatitude() < b.getLatitude())) ||
                ((b.getLatitude() <= point.getLatitude()) && (point.getLatitude() < a.getLatitude()))) {
                double vt = (point.getLatitude() - a.getLatitude()) / (b.getLatitude() - a.getLatitude());
                if (point.getLongitude() < a.getLongitude() + vt * (b.getLongitude() - a.getLongitude())) {
                    crossings++;
                }
            }
        }
        return (crossings % 2 == 1);
    }

private:
    std::vector<GeoCoordinateStruct> m_points;  // 多边形的点
};

//定义散列点类
class HashPointStruct {
public:
	//默认构造函数
	HashPointStruct() = default;

	//构造函数,默认参考点为(0,0)
    HashPointStruct(std::vector<GeoCoordinateStruct> coordinates, GeoCoordinateStruct reference = GeoCoordinateStruct(0, 0)) {
		m_coordinates = coordinates;
		m_count = coordinates.size();
		m_reference = reference;

        //计算散列点中心点坐标
        m_center = HashPointStruct::calculateCenterPoint(m_coordinates);
	}

    //返回散列点个数
    int getCount() { return m_count; }

    //返回散列点中心点坐标
    GeoCoordinateStruct getCenter() { return m_center; }

    //计算散列点中心点坐标
    static GeoCoordinateStruct calculateCenterPoint(std::vector<GeoCoordinateStruct> coordinates) {
		//计算散列点中心点坐标
		double longitude = 0;
		double latitude = 0;
        for (auto c : coordinates) {
			longitude += c.getLongitude();
			latitude += c.getLatitude();
		}
		longitude /= coordinates.size();
		latitude /= coordinates.size();

		return GeoCoordinateStruct(longitude, latitude);
	}

    //计算矩形区域
 //   RectangleStruct calculateRectangle() {
 //       //计算直线
 //       LineStruct line1(m_center, m_reference);

 //       //计算方位角
 //       double bearing = m_reference.calculateBearing(m_center);



	//	return rectangle;
	//}

private:
    std::vector<GeoCoordinateStruct> m_coordinates;	//散列坐标列表
	GeoCoordinateStruct m_center;	//散列坐标中心点坐标
    int m_count;            //散列点个数
    GeoCoordinateStruct m_reference;    //参考点
};

//功能：计算弹药投放矩形区域
//输入1：一个多边形阵地区域
//输入2：敌方若干实体坐标
//输入3：我方若干实体坐标
//输出：弹药投放矩形区域
RectangleStruct calculateRectArea(std::vector<GeoCoordinateStruct> polygon, std::vector<GeoCoordinateStruct> enemy, std::vector<GeoCoordinateStruct> ours) {
    //1.判断敌方实体是否在阵地区域内，并计算敌方实体中心点坐标
    std::vector<GeoCoordinateStruct> enemyInPolygon;
    for (auto e : enemy) {
        if (PolygonStruct(polygon).contains(e)) {
            enemyInPolygon.push_back(e);
        }
    }

    //如果enemyInPolygon为空，让enemyInPolygon等于polygon
    if (enemyInPolygon.empty()) {
        enemyInPolygon = polygon;
    }

    GeoCoordinateStruct enemyCenter = calculateCoordsCenter(enemyInPolygon);

    //2.计算我方所有实体的中心点坐标
    GeoCoordinateStruct oursCenter = calculateCoordsCenter(ours);

    //3.以我方中心点和敌方中心点为两个点确定一条线，计算出斜率和垂直于该线的斜率
    LineStruct line(oursCenter, enemyCenter);
    double slope = line.getSlope();
    double perpendicularSlope = line.calculatePerpendicularSlope();

    //4.遍历所有敌方实体，分别计算出过该点两条直线的截距并分别保存截距最大和最小的两条直线
    LineStruct maxSlopeLine(slope, enemyInPolygon[0]);
    LineStruct minSlopeLine(slope, enemyInPolygon[0]);
    LineStruct maxPerpendicularSlopeLine(perpendicularSlope, enemyInPolygon[0]);
    LineStruct minPerpendicularSlopeLine(perpendicularSlope, enemyInPolygon[0]);

    for (auto e : enemyInPolygon) {
        LineStruct l1(slope, e);
        if (l1.getYIntercept() > maxSlopeLine.getYIntercept()) {
            maxSlopeLine = l1;
        }

        if (l1.getYIntercept() < minSlopeLine.getYIntercept()) {
            minSlopeLine = l1;
        }

        LineStruct l2(perpendicularSlope, e);
        if (l2.getYIntercept() > maxPerpendicularSlopeLine.getYIntercept()) {
            maxPerpendicularSlopeLine = l2;
        }

        if (l2.getYIntercept() < minPerpendicularSlopeLine.getYIntercept()) {
            minPerpendicularSlopeLine = l2;
        }
    }

    //5.根据4条直线确定矩形区域
    std::vector<LineStruct> lines = { maxSlopeLine, maxPerpendicularSlopeLine, minSlopeLine, minPerpendicularSlopeLine };
    RectangleStruct rectangle(lines, oursCenter);

    return rectangle;
}

//功能：【榴弹炮火力分配算法】根据射击模式不同对阵地矩形区域进行分割，计算弹药投放落点
//输入1：阵地矩形区域
//输入2：射击模式(1,连单独射击；2,营重叠射击)
//输出：弹药投放落点
std::vector<GeoCoordinateStruct> calculateDropPoints(RectangleStruct rectangle, int mode) {
    //1.矩形正面和纵深以50米为单位向上归整
    double front = ceil(rectangle.getFront() / 50) * 50;
    double depth = ceil(rectangle.getDepth() / 50) * 50;

    //2.根据射击模式计算弹药投放落点
	std::vector<GeoCoordinateStruct> dropPoints;
    if (mode == 1) {
        //连单独射击
        if (front <= 50 && depth <= 200) {
            //以矩形中心点为落点
            dropPoints.push_back(rectangle.calculateCenter());
        }
        else if (front > 50 && front <= 100 && depth <= 200) {
            //分割矩形计算落点
            std::vector<RectangleStruct> rectangles = rectangle.split(6, 1);
            for (auto r : rectangles) {
                dropPoints.push_back(r.calculateCenter());
            }
        }
        else if (front > 100 || depth > 200) {
            //矩形分段计算
            int m = ceil(front / 100);
            int n = ceil(depth / 200);
            std::vector<RectangleStruct> rectangles = rectangle.split(m, n);

            //打印分割后的矩形
            for (auto r : rectangles) {
                std::cout << "矩形区域的四个顶点坐标：" << std::endl;
				std::cout << r.getTopLeft().toDMSString() << std::endl;
				std::cout << r.getTopRight().toDMSString() << std::endl;
				std::cout << r.getBottomRight().toDMSString() << std::endl;
				std::cout << r.getBottomLeft().toDMSString() << std::endl;
			}

            for (auto r : rectangles) {
                //分割矩形计算落点
                std::vector<RectangleStruct> subRectangles = r.split(6, 1);
                for (auto sr : subRectangles) {
                    dropPoints.push_back(sr.calculateCenter());
                }
            }
        }
        else {
            ;
        }
	}
    else if (mode == 2) {
        //营重叠射击
        if (front <= 100 && depth <= 300) {
            //以矩形中心点为落点
            dropPoints.push_back(rectangle.calculateCenter());
        }
        else if (front > 100 && front <= 300 && depth <= 300) {
            //分割矩形计算落点
            std::vector<RectangleStruct> rectangles = rectangle.split(6, 1);
            for (auto r : rectangles) {
                dropPoints.push_back(r.calculateCenter());
            }
        }
        else if (front > 100 && front <= 300 && depth > 300 && depth <= 600) {
            //分割矩形计算落点
            std::vector<RectangleStruct> rectangles = rectangle.split(6, 3);
            for (auto r : rectangles) {
                dropPoints.push_back(r.calculateCenter());
            }
        }
        else if (front > 300 || depth > 600) {
			//矩形分段计算
			int m = ceil(front / 300);
			int n = ceil(depth / 600);
			std::vector<RectangleStruct> rectangles = rectangle.split(m, n);
            for (auto r : rectangles) {
				//分割矩形计算落点
				std::vector<RectangleStruct> subRectangles = r.split(6, 3);
                for (auto sr : subRectangles) {
					dropPoints.push_back(sr.calculateCenter());
				}
			}
		}
        else {
			;
		}
	}
	return dropPoints;
}

//功能：【远程火箭炮火力分配算法】
std::vector<GeoCoordinateStruct> calculateDropPoints(RectangleStruct rectangle)
{
    //1.矩形正面和纵深以50米为单位向上归整
    double front = ceil(rectangle.getFront() / 50) * 50;
    double depth = ceil(rectangle.getDepth() / 50) * 50;

    //2.计算弹药投放落点
    std::vector<GeoCoordinateStruct> dropPoints;
    if (front <= 150 && depth <= 150) {
        //以矩形中心点为落点
        dropPoints.push_back(rectangle.calculateCenter());
    }
    else if (front > 150 || depth > 150) {
        //矩形分段计算
        int m = ceil(front / 150);
        int n = ceil(depth / 150);
        std::vector<RectangleStruct> rectangles = rectangle.split(m, n);

        //打印分割后的矩形
        for (auto r : rectangles) {
            std::cout << "矩形区域的四个顶点坐标：" << std::endl;
            std::cout << r.getTopLeft().toDMSString() << std::endl;
            std::cout << r.getTopRight().toDMSString() << std::endl;
            std::cout << r.getBottomRight().toDMSString() << std::endl;
            std::cout << r.getBottomLeft().toDMSString() << std::endl;
        }

        for (auto r : rectangles) {
            dropPoints.push_back(r.calculateCenter());
        }
    }
    else {
        ;
    }

    return dropPoints;
}

//测试【榴弹炮火力分配算法】
void TestCalculateDropPoints1()
{
    //1.输入阵地区域
    std::vector<GeoCoordinateStruct> polygon;
    polygon.push_back(GeoCoordinateStruct(91.789815, 28.175279));
    polygon.push_back(GeoCoordinateStruct(91.792692, 28.125399));
    polygon.push_back(GeoCoordinateStruct(91.841881, 28.112551));
    polygon.push_back(GeoCoordinateStruct(91.86123, 28.170731));

    //2.输入敌方实体坐标
    std::vector<GeoCoordinateStruct> enemy;
    enemy.push_back(GeoCoordinateStruct(91.812, 28.158));
    enemy.push_back(GeoCoordinateStruct(91.822, 28.158));
    enemy.push_back(GeoCoordinateStruct(91.831, 28.159));
    enemy.push_back(GeoCoordinateStruct(91.811, 28.144));
    enemy.push_back(GeoCoordinateStruct(91.821, 28.143));
    enemy.push_back(GeoCoordinateStruct(91.831, 28.142));

    //3.输入我方实体坐标
    std::vector<GeoCoordinateStruct> ours;
    ours.push_back(GeoCoordinateStruct(91.967, 27.936));
    ours.push_back(GeoCoordinateStruct(91.974, 27.936));
    ours.push_back(GeoCoordinateStruct(91.980, 27.937));
    ours.push_back(GeoCoordinateStruct(91.987, 27.937));
    ours.push_back(GeoCoordinateStruct(91.994, 27.938));
    ours.push_back(GeoCoordinateStruct(92.002, 27.938));

    //4.计算弹药投放矩形区域，并采用度分秒格式输出矩形的四个顶点
    RectangleStruct rectangle = calculateRectArea(polygon, enemy, ours);
    std::cout << "矩形区域的四个顶点坐标：" << std::endl;
    std::cout << rectangle.getTopLeft().toDMSString() << std::endl;
    std::cout << rectangle.getTopRight().toDMSString() << std::endl;
    std::cout << rectangle.getBottomRight().toDMSString() << std::endl;
    std::cout << rectangle.getBottomLeft().toDMSString() << std::endl;
    //打印正面和纵深
    std::cout << "矩形区域的正面和纵深：" << std::endl;
    std::cout << "正面：" << rectangle.getFront() << "米" << std::endl;
    std::cout << "纵深：" << rectangle.getDepth() << "米" << std::endl;

    //5.计算弹药投放落点
    std::vector<GeoCoordinateStruct> dropPoints = calculateDropPoints(rectangle, 1);
    std::cout << "弹药投放落点个数：" << dropPoints.size() << std::endl;
    std::cout << "弹药投放落点坐标：" << std::endl;
    for (auto point : dropPoints) {
        std::cout << point.toDMSString() << std::endl;
    }
}

//测试【远程火箭炮火力分配算法】
void TestCalculateDropPoints2()
{
    //1.输入阵地区域
    std::vector<GeoCoordinateStruct> polygon;
    polygon.push_back(GeoCoordinateStruct(91.789815, 28.175279));
    polygon.push_back(GeoCoordinateStruct(91.792692, 28.125399));
    polygon.push_back(GeoCoordinateStruct(91.841881, 28.112551));
    polygon.push_back(GeoCoordinateStruct(91.86123, 28.170731));

    //2.输入敌方实体坐标
    std::vector<GeoCoordinateStruct> enemy;
    enemy.push_back(GeoCoordinateStruct(91.812, 28.158));
    enemy.push_back(GeoCoordinateStruct(91.822, 28.158));
    enemy.push_back(GeoCoordinateStruct(91.831, 28.159));
    enemy.push_back(GeoCoordinateStruct(91.811, 28.144));
    enemy.push_back(GeoCoordinateStruct(91.821, 28.143));
    enemy.push_back(GeoCoordinateStruct(91.831, 28.142));

    //3.输入我方实体坐标
    std::vector<GeoCoordinateStruct> ours;
    ours.push_back(GeoCoordinateStruct(91.967, 27.936));
    ours.push_back(GeoCoordinateStruct(91.974, 27.936));
    ours.push_back(GeoCoordinateStruct(91.980, 27.937));
    ours.push_back(GeoCoordinateStruct(91.987, 27.937));
    ours.push_back(GeoCoordinateStruct(91.994, 27.938));
    ours.push_back(GeoCoordinateStruct(92.002, 27.938));

    //4.计算弹药投放矩形区域，并采用度分秒格式输出矩形的四个顶点
    RectangleStruct rectangle = calculateRectArea(polygon, enemy, ours);
    std::cout << "矩形区域的四个顶点坐标：" << std::endl;
    std::cout << rectangle.getTopLeft().toDMSString() << std::endl;
    std::cout << rectangle.getTopRight().toDMSString() << std::endl;
    std::cout << rectangle.getBottomRight().toDMSString() << std::endl;
    std::cout << rectangle.getBottomLeft().toDMSString() << std::endl;
    //打印正面和纵深
    std::cout << "矩形区域的正面和纵深：" << std::endl;
    std::cout << "正面：" << rectangle.getFront() << "米" << std::endl;
    std::cout << "纵深：" << rectangle.getDepth() << "米" << std::endl;

    //5.计算弹药投放落点
    std::vector<GeoCoordinateStruct> dropPoints = calculateDropPoints(rectangle);
    std::cout << "弹药投放落点个数：" << dropPoints.size() << std::endl;
    std::cout << "弹药投放落点坐标：" << std::endl;
    for (auto point : dropPoints) {
        std::cout << point.toDMSString() << std::endl;
    }
}

//测试构造矩形
void TestConstructRectangle()
{
    GeoCoordinateStruct center(91.969456, 27.176189);
    double bearing = 35;
    double front = 2000;
    double depth = 1000;

    RectangleStruct rectangle(center, bearing, front, depth);
    
    //打印矩形
    std::cout << "矩形区域的四个顶点坐标：" << std::endl;
    std::cout << rectangle.getTopLeft().toDMSString() << std::endl;
    std::cout << rectangle.getTopRight().toDMSString() << std::endl;
    std::cout << rectangle.getBottomRight().toDMSString() << std::endl;
    std::cout << rectangle.getBottomLeft().toDMSString() << std::endl;

    //分割成2*2
    std::vector<RectangleStruct> rectangles = rectangle.split(2, 2);

    //打印分割后的矩形中心点
    std::cout << "分割后的矩形中心点坐标：" << std::endl;
    for (auto r : rectangles) {
	    std::cout << r.calculateCenter().toDMSString() << std::endl;
    }
}
